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SUMMARY 

A number of approaches for improving the accuracy of incompressible, 
steady-state flow calculations are examined. Two Improved differencing 
schemes, Quadratic Upstream Interpolation for Convective Kinematics (QUICK) and 
Skew-Upwind Differencing (SUD), are applied to the convective terms In the 
Navler-Stokes equations and compared with results obtained using hybrid differ- 
encing. In a number of test calculations, It Is Illustrated that no single 
scheme exhibits superior performance for all flow situations. However, both 
SUD and QUICK are shown to be generally more accurate than hybrid differencing. 


INTRODUCTION 

Three-dimensional calculations of combustor flow fields are currently 
imposing severe demands on the computer hardware capabilities and computing 
budgets of gas turbine manufacturers. One of the main reasons for this relates 
to the large number of complex physical processes occurring in the combustor. 
Airflow, fuel spray, reaction kinetics, flame radiation, and turbulence must 
be modeled and the related differential equations solved. Obtaining accurate 
solutions to these modeled equations entails overcoming another difficulty. 
Current combustor codes are, generally, based on the SIMPLE algorithm developed 
by Patankar and Spalding (ref. 1). Frequently used In conjunction with this 
solution algorithm Is hybrid differencing to approximate the convective terms 
In the governing equations. In most practical calculations this results In 
the use of first order accurate upwind differencing for most of the flow field. 
The overly dissipative solutions obtained can mask Important flow field fea- 
tures (ref. 2). Ideally, this can be overcome through the use of additional 
grid points. In three dimensional calculations, however, this soon becomes 
Impractical due to the large number of physical processes that must be modeled. 

To alleviate this problem, NASA has Initiated a program to Identify and 
Incorporate an Improved accuracy differencing scheme Into a combustor perform- 
ance code. Under a portion of this program a competitive contract was awarded 
to United Technologies (with A.D. Gosman as consultant) to examine and Imple- 
ment a variety of Improved accuracy differencing schemes. The schemes examined 
Include Quadratic Upstream Interpolation for Convective Kinematics (QUICK) and 
variants of Skewed-Upwlnd Differencing (SUD). The variants of SUD Involve dif- 
ferent procedures to ensure that nonphysical oscillations do not occur In the 
solution. 




The purpose of this paper is to draw some general conclusions concerning 
the accuracy of the convective differencing schemes studied. A number of non- 
turbulent test calculations are made to illustrate the accuracy and stability 
(convergence) characteristics of the various schemes. In general, these test 
calculations indicate that both QUICK and SUDS are more accurate than hybrid 
differencing although these schemes greatly slow convergence of the governing 
equations . 


MATHEMATICAL FORMULATION 


The numerical code employed in this study provides the capability to 
analyze steady-state, two-dimensional, elliptic, turbulent flows. Only the 
pertinent features of the code shall be reviewed here as further details are 
available in the open literature. 


The governing equations, written in tensor notation, include: 
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The numerical code used in these studies can solve the two additional 
equations relating to the turbulence model of reference 3. However, the test 
calculations are meant to focus on the accuracy with which the momentum equa- 
tions can be solved, therefore, the viscosity was fixed at a laminar value. 

The governing equations are discretized through integral analysis or the 
"finite volume" method following the procedure of Gosman and Lai, (ref. 4). 
Integrals across each computational cell face are evaluated using the Mean 
Value Theorem on the grid illustrated in figure 1(a). Approximations for the 
convective and diffusive terms are then made depending on the differencing 
procedure used. The pressure gradient and diffusive terms are central differ- 
enced while the differencing of the convective terms will be noted in the next 
section. The resulting fluxes across each cell face are then summed. This 
equation when arranged into a substitution formula for the variable being 
solved (for example 4 > p is the variable being solved for at point P) 
becomes : 


a p 4> p =^a i 4 * i + Su (3) 

K K W 1 

r 

where denotes summation of the neighbors of P. 


2 


CONVECTIVE DIFFERENCING 


Hybrid Differencing 

The practice most commonly used, at present, is to employ hybrid differ- 
encing to approximate the convective terms. This involves the use of second 
order accurate central differencing when the cell Peclet number (Pe c ) is less 
than an absolute value of 2, while first order accurate upwind differencing is 
used when |Pe c | > 2. The major advantage of this scheme is the uncondi- 
tionally bounded or oscillation free solutions it provides. 

For example, if solving for <t> on the staggered grid (fig. 1(a)), one of 
the convective terms becomes: 


„ (u<|>) - ( Ud>) 

aud) _ 'e v w 

ax " ax 

where for central differencing (uniform grid spacing) - 

(u+) e = 2 u e (<t, E + V 

<“♦>« " 2 V+W * V 


for upwind differencing (u > 0) - 

( u ^) e = Vp 

(u ^w = 2 U w<*W + V 


(4) 


QUICK Differencing 

This scheme, developed by Leonard (ref. 5), is second order accurate and 
not unconditionally bounded. (The scheme may produce nonphysical oscillations 
in the solution.) It approximates convective terms using an upstream biased 
quadratic interpolation. For example, using the grid in figure 1(b) and 
approximating equation (4) as before: 


For u > 0 


U e 

< u *>w= U w 


8 *W * 4 *P + 8 +El 


" 8 *WW + 4 *W + 8 *P 


Skew Upwind Differencing (SUD) 

This scheme, developed by Raithby (ref. 6), Is first order accurate and, 
as with QUICK, is not unconditionally bounded. While SUD is formally the same 
order accuracy as upwind differencing, its truncation error is less than upwind. 
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SUO attains this higher accuracy fay differencing in an upwind manner along the 
flow streamlines. Each streamline is defined by the velocity vector at each 
grid boundary (fig. 1(c)). The upstream value of the variable to be calculated 
is then obtained by a back projection of the velocity vector and simply inter- 
polating between the two neighboring values. For example, using the grid of 
figure 1(c) and solving equation (4) as in the previous two examples: 

For u > 0 and v > 0 


Vp 


0*>w - u w o - «> * w + v* 


sw 


where 


a = minimum of (1, v/2u) 

Flux Blending 

The concept of flux-blending is analogous to the "Flux Corrected Trans- 
port" (FCT) technique of Boris and Books (ref. 7). The procedures employed 
here were developed by Gosman, Lai, and Peric and are detailed in reference 8. 
In general, the flux blended schemes employ a weighted mean of a bounded (but 
low order accuracy) differencing scheme and an unbounded, more accurate scheme. 
The main factor is to blend as little of the lesser accurate scheme as possible 
while still maintaining a properly bounded solution. The two differencing 
schemes blended involve upwind differencing and the more accurate SUD. 


BSUDS1 

This procedure, "Bounded Skew Upwind Differencing Scheme 1" blends upwind 
and SUD in proportions sufficient to suppress negative coefficients in 
equation (3). This is a sufficient condition for a bounded solution, however, 
it tends to be more dissipative than necessary. 


BSUDS2 

This procedure "Bounded Skew Upwind Differencing Scheme 2" blends upwind 
and SUD in proportions ensuring that when negative coefficients occur, their 
contribution is below the level that would cause the solution to be physically 
unrealistic. This procedure is Iterative and starts from an initial, totally 
skew differenced estimate. If the calculated variable has a value no greater 
or lesser than that of its neighbors, then the solution is bounded and no 
blending is performed. If the solution is out of the range of neighboring 
values, then blending is performed. In the extreme, this blending would result 
In upwind differencing. The use of neighboring values as limits in determin- 
ing the "boundedness" of the solution is only valid when the equation being 
solved lacks source terms. However, the momentum equations contain signifi- 
cant source contributions. The Implications of this are still being studied, 
but the results reported in later sections shall demonstrate that the use of 
neighboring values as physical limits provides highly accurate results. 
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This bounding procedure, while simple In concept, Is difficult to apply to 
an Iterative solution scheme. If an Initial SUD calculation was made and then 
the coefficients were updated for bounding and the equation solved a second 
time, the computational time required for one Iteration would be approximately 
doubled. To reduce this computational overhead BSUDS2 calculations were typi- 
cally restarted from a BSUDS1 calculation with the bounding evaluated based on 
the previous Iteration values. This results In some "unboundedness" when the 
equations are not fully converged, however, the final result should be bounded. 


Solution Algorithm 

As solution algorithm for solving the governing equations will be only 
briefly reviewed here (refs. 4 and 9) are strongly recommended for further 
detal 1 s . 

Once the momentum equations are approximated on the staggered grid, these 
equations must be solved In a process Insuring that all governing equations are 
satisfied. In the SIMPLE algorithm, each momentum equation Is sequentially 
solved using a guessed or old pressure field from the previous Iteration. A 
pressure correction equation Is then solved and the values of velocity and 
pressure are revised to more closely satisfy continuity. Iteration on this 
procedure Is then continued until all equations are satisfied to a low normal- 
ized residual level (typically 10 -2 or 10~ 4 for the calculations herein 
reported) . 


Scalar Transport Test Calculation 

Following the procedure detailed In reference 8, a preliminary assessment 
of the accuracy associated with the various differencing schemes was made solv- 
ing the following scalar transport equation: 

iJJ£ + _ Q 

3x ay 


for 

U = constant 

V = constant = U tan 0 

where 

0 = flow angle 

A normal distribution of $ was Imposed along a streamline centered 
coordinate (fig. 2). This provided the boundary conditions for <j> and a 
single point (<t>p) calculation could then be made using the differencing schemes 
described earlier. 
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RESULTS AND DISCUSSION 


The accuracy of the previously described convective differencing schemes 
Is Illustrated through a series of test calculations. The first test problem 
Is a single cell calculation of scalar transport. Following this example a 
variety of laminar flow calculations are made to examine whether the conclu- 
sions that were drawn from the scalar calculation remain valid for the full 
momentum equations. The effect that each of these schemes has on convergence 
Is examined in the final section. 


Scalar Transport Test Calculation 

A straight-forward approach to study the accuracy of various differencing 
schemes Is to solve an equation describing the transport of a scalar with no 
diffusion or source terms. This problem has been described In the previous 
section and Is graphically Illustrated (fig. 2). The exact answer to this 
problem Is that the scalar, <j>, remains equal to 1.0 along the streamline 
Intersecting point P. 

Solutions for the scalar transport equation as a function of flow angle 
for the various differencing schemes are displayed in figure 3. All schemes 
agree with the exact solution at a flow angle of zero, but departing from this 
each scheme displays some degree of error relating to numerical smearing or 
diffusion. Upwind differencing displays the greatest amount of numerical dif- 
fusion with the largest errors occurring at 45°. QUICK displays a similar 
behavior but the overall level of error Is much less. SUD displays an error 
maximum around 15° but tends to zero at angles approaching 45°. BSUDS1 dis- 
plays more error than SUD with a maximum error around 22°. All the first order 
schemes, upwind, SUD, and BSUDS1 , display similar error levels for flow angles 
less than about 5° although these errors are generally small. QUICK displays 
superior performance to the best of the first order schemes, SUD, at flow 
angles up to about 15°. Above a flow angle of 15° SUD Is generally more accu- 
rate than any of the other schemes. 

The penalty to be paid for ensuring against nonphysical oscillations using 
the first flux blending scheme can be seen by comparing the SUD and BSUDS1 
results In figure 3. (The first flux blending scheme requires that all advec- 
tlon coefficients (eq. (3)) be positive as detailed In the earlier section.) 

For example, at 30° BSUDS1 displays about 6 percent error while SUD displays 
less than one percent. The BSUDS2 scheme Is not readily suitable to test this 
problem, but logically one can expect this scheme to exhibit an accuracy 
between SUD and BSUDS1 . 

It Is Important to note that these test calculations were made on a uni- 
formly spaced grid. The accuracy of the SUD calculations would change as the 
cell aspect ratio (AX/AY) was varied. Figure 4 displays an example of how cell 
aspect ratio affects the accuracy of SUD. (This figure Is shown for Illustra- 
tion purposes only - the Indicated accuracy Is dependent on the formulation of 
the scalar calculation.) For Increasing aspect ratios above 1, SUD Is more 
accurate at flow angles only slightly skewed to the grid. For example, at a 
cell aspect ratio of 2, SUD has an optimal performance around a flow angle of 
26°. At an aspect ratio of 1, the best performance of SUD Is around 45°. 


Due to the impact of cell aspect ratio on SUD accuracy all of the test 
calculations were made on a uniform mesh or the mesh was uniform in the impor- 
tant flow regions and allowed to be nonuniform only far downstream. 


Laminar Flow Test Calculations 

To examine whether the scalar calculations remain a valid test of convec- 
tive accuracy when the full momentum equations are solved, a series of laminar 
flow calculations were made. Laminar test calculations were chosen over tur- 
bulent flow calculations due to the fact that numerical inaccuracy should 
exhibit itself as numerical smearing. The diffusion added by a turbulence 
model would obsecure the differences In convective accuracy. In laminar flow 
calculations the more accurate solutions will display steeper velocity 
gradients. 


Driven Cavity Test Calculations 

The first laminar flow calculation is of the driven cavity flow field. 

This geometry Is frequently used to assess numerical accuracy (refs. 10 to 13). 
Fine grid calculations (refs. 11 and 12) are available as a "standard" for 
accuracy comparisons. Figure 5 displays a typical flow field calculation for 
this geometry. The moving upper wall (not shown) imposes a circulation in the 
cavity. The calculations were made for a Reynold's number of 1000 based on 
cavity height. The grid points in these calculations were uniformly spaced. 

Calculated velocity profiles through the center of the circulation vortex 
are shown In figure 6. In this figure the results of progressively finer grids 
are compared against the "standard". The hybrid results are slowly tending 
toward the more accurate results and yet, even the 80 by 80 calculations 
exhibit numerical inaccuracy. 

The velocity peak seen in figure 6 around Y/L = 0.2 is used as a measure 
of accuracy in figure 7. Although it is difficult to quantify the accuracy of 
an entire two-dimensional flow field by a single point, this point does seem 
representative of the accuracy. In figure 7 the velocity peak is shown versus 
the number of grid points used in the hybrid, QUICK, and BSUDS2 calculations. 
The line at U/Uw = 0.4 indicates the "standard" result. The QUICK results 
indicate a strong response to grid refinement such that at 6400 grid points the 
QUICK results match the "standard". The first order schemes hybrid and BSUDS2 
appear to be highly Inaccurate with the 6400 grid point calculations being only 
as accurate as the 400 grid point calculations of QUICK. 

The poor performance of the BSUDS2 calculations Is surprising in view of 
the scalar calculations previously shown. Calculations were made using 
unbounded SUD and were found to be indistinguishable from the BSUDS2 calcula- 
tions. Therefo.e, the flux blended procedure can not account for the poor 
performance. 


Laminar Flow Calculations With Various Inlet Flow Angles 

A second series of laminar flow calculations were made to further explore 
the flow angle/accuracy dependance. The geometry employed in these calculat- 
ions is shown in figure 8. The inlet flow angle was varied from 40 to 0°. 
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All calculations were for a Reynold's number of 1000 based on Inlet height. 
Grid points were arranged with uniform spacing In the forward portion of the 
flow field (aspect ratio equal to 1) while the mesh was nonuniform far 
downstream. 


Inlet Flow Angle of 40° 

Calculations made with two different grid meshes are shown for a X/H of 
0.5 In figure 9. The coarse mesh calculations (fig. 9(a)) show that both QUICK 
and BSUDS2 provide steeper gradients than the hybrid calculation. The BSUDS2 
results appear to be overall more accurate especially when these results are 
compared to the fine mesh calculations (fig. 9(b)). The fine mesh calculations 
(fig. 9(b)) again display the poor performance of hybrid differencing and the 
higher accuracy of QUICK and BSUDS2. In contrast to the coarse mesh results 
qUICK and BSU0S2 appear to be of nearly similar accuracy. 

The performance of BSU0S1 and BSUDS2 for the fine mesh calculations can be 
seen In figure 9(c). The results for these schemes are quite similar with the 
BSUDS2 exhibiting slightly steeper velocity gradients, and thereby slightly 
greater accuracy. On a coarse mesh (not shown) the results were similar. 


Inlet Flow Angle of 25° 

Velocity profiles at X/H = 0.5 for an Inlet flow angle of 25° are dis- 
played In figure 10. The coarse grid results (fig. 10(a)) Indicates that 
hybrid differencing again smears the velocity profile more than the QUICK or 
BSUDS2 results. The fine mesh results (fig. 10(b)) Indicate smaller differ- 
ences between all three schemes. The hybrid calculation smears the velocity 
profile more than the other two calculations, but the Indicated error Is less 
than that in the 40° test calculations of figure 9. The QUICK results exhibit 
a nonphysical oscillation In the velocity profile around Y/H = 0.5 while the 
BSU0S2 result displays a smooth but accurate velocity profile. 

The comparison between the BSU0S1 and the BSUDS2 calculations Is displayed 
In figure 10(c). The results are similar to the 40° comparison (fig. 9(c)) 
with BSUDS2 displaying only slightly Improved accuracy. This result Is sur- 
prising In light of the scalar calculations of figure 3. In those calculations 
there was a large variation In accuracy between BSUDS1 and BSUDS2 at the two 
different flow angles. This will be more fully discussed In the summary 
section. 


Inlet Flow Angle of 0° 

Velocity profiles at X/H = 1.0 for an Inlet flow angle of 0° are dis- 
played In figure 11. Both the coarse mesh results (fig. 11(a)) and the fine 
mesh results (fig. 11(b)) show very little effect of the convective operator. 
This Is expected since the scalar transport test calculations Indicated that 
all schemes are essentially exactly correct for a 0° flow angle. 
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Summary of Accuracy Tests 

From the series of flow calculations made, no single differencing scheme 
displayed superior accurac / In all calculations. The scalar calculation demon- 
strated that performance , the differencing scheme was related to flow angle. 
The laminar flow calculations Indicated that other factors must be equally 
Important. For Instance, a possible explanation of the relatively poor per- 
formance of BSUDS2 In the driven cavity calculation versus the other flow cal- 
culations can be advanced from consideration of the gradients In the calculated 
variable. In the Inclined Inlet flow calculations, the linear variation of 
velocity Issuing from the Inlet can be well represented by any form of SUD. 

The normal distribution used In the scalar calculations can also be fairly 
accurately accommodated by SUD. The driven cavity exhibits velocity gradients 
which are dissimilar to either of the previous examples. The moving wall 
establishes a steep velocity gradient which can, apparently, only be fairly 
well represented by QUICK. However this driven cavity calculation may be some- 
what unrepresentative of combustor-type geometries. The Inclined Inlet flow 
field has features that are highly similar to combustors, suggesting that SUD 
can Le used In combustor calculations to provide some Improved accuracy. 


Convergence Tests 

Computational times to converge the system of governing equations Is shown 
In table I for three different convective schemes. These calculations were all 
made on a Cray Is using the SIMPLE algorithm for the driven cavity flow geome- 
try. In general, either BSUDS2 or QUICK required more time to converge than 
hybrid. The slowest performer was QUICK with times In the range of three to 
five times slower to converge than hybrid differencing. BSUDS2 was usually 
only slightly slower than hybrid. The other variants of SUD gave results 
comparable to the BSUDS2 timings. It should be noted that the choice of under- 
relaxation has a strong Influence on the convergence rate and that the highest 
under-relaxation that would promote stabile convergence was used. However, 
extensive optimization was not done and It may be possible to achieve slightly 
faster or slower convergence times. 


CONCLUSIONS 

1. The accuracy of BSUDS2 and QUICK was generally greater than hybrid 
differencing for the series of calculations reported here. In the driven 
cavity series of calculations BSUDS2 was about as Inaccurate as hybrid dif- 
ferencing while QUICK displayed a high level of accuracy. In a series of 
calculations more representative of combustor-type geometries (Inclined 
Inlet), BSUDS2 and QUICK were highly accurate while hybrid differencing was 
very Inaccurate. 

2. The bounding procedure used In BSUDS2 produced solutions frr-e of 
nonphysical oscillations while maintaining a high level of accuracy. 

3. Both QUICK and SUD generally required more time to converge than the 
hybrid calculations. This became more pronounced In fine mesh calculations. 
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TABLE I. - CONVERGENCE TIMES FOR 
THREE CONVECTIVE DIFFERENCING 
SCHEMES USING THE SIMPLE 
ALGORITHM. DRIVEN CAVITY 
FLOW CALCULATIONS 


Number 

of 

nodes 

Hybrid 

BSUDS2 

QUICK 

CPU 

time, sec 

100 

6.1 

11.9 

7.6 

400 

17.5 

18.0 

47.1 

1600 

118.5 

118.2 

524.7 

6400 

1019.2 

1517.6 

5551.0 
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Figure 2. - Illustration of the scalar transport test calculation* 
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Figure 3* - Results of the scalar transport test 
calculations shown versus flow angle. 






CELL ASPECT RATIO, AX/AY 

Figure 4. - Optimum flow angle for skew differ- 
encing (SUDI shown for various computational 
cell aspect ratios. Shaded region Indicates 
approximately less than one percent Inaccuracy 
In solution of scalar transport test calculation. 



Figure 5. - Velocity vectors for a typical driven cavity calculation. 
The upper moving wall Is removed from the plot Flow Re * 1000, 



Figure 6. - Velocity profiles through the center of the 
circulation vortex for the driven cavity calculations. 
Shown are the results of hybrid calculations compared 
with high accuracy "standard" results. 
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Figure 7. - Velocity peak value plotted against number of grid points for the 
driven cavity flow calculations. 
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Figure 8, - Geometry Illustration for laminar flow calculations 
employing various Inlet flow angles. 








Figure la - Laminar flew/ test calculations 
comparing hybrid, QUICK and BSU0S2. 
Inlet flow angle of 25°; geometry of 
figure 8. X/H-a5, 
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(c) Fine grid (58x38) calculations compar- 
ing BSUDS1 and BSUDS* 

Figure la - Concluded. 






